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Abstract 

Motivated by a controversy over the correct results derived from the dynamic 

renormalization group (DRG) analysis of the non linear molecular beam epitaxy 

(MBE) equation, a self-consistent expansion (SCE) for the non linear MBE theory is 

considered. The scaling exponents are obtained for spatially correlated noise of the 

general form D(r -f ',t-t')= 2D If -r 1 S(t-t'). I find a lower critical 

dimension d c (p) = 4 + 2p , above which the linear MBE solution appears. Below the 

lower critical dimension a p-dependent strong-coupling solution is found. These 
results help to resolve the controversy over the correct exponents that describe non 
linear MBE, using a reliable method that proved itself in the past by giving reasonable 
results for the strong coupling regime of the Kardar-Parisi- Zhang (KPZ) system (for 
d > 1 ), where DRG failed to do so. 



The field of disorderly surface growth has received much attention during the 
last two decades. Special effort has been focused on relating discrete microscopic 
growth models with their corresponding continuum field theories [1]. The first 
continuum equation used to study the growth of interfaces by particle deposition was 
the Edwards- Wilkinson model (EW) [2] that describes the dynamics of the interface 
by a noise driven diffusion equation. This model actually describes the microscopic 
process known as random deposition (RD) with surface relaxation, and together they 
form a distinct universality class in growth phenomena. However, an extension to this 
model was needed because of the nonlinear character of many deposition processes, 
such as ballistic deposition (BD), solid-on- solid deposition (SOS) and Eden growth. 
The first extension of the EW equation to include nonlinear terms was proposed by 
Kardar, Parisi and Zhang (KPZ) [3], who suggested the addition of a nonlinear term 
proportional to the square of the height gradient. The success of the KPZ equation in 
describing deposition phenomena motivated many researchers to develop a continuum 
growth model relevant for the technologically important molecular beam epitaxy 
(MBE) process [4-10]. The physical mechanism that distinguishes MBE from 
previously discussed growth processes is the surface diffusion of the deposited 
particles. It is well known that in the temperature range of MBE growth, desorption of 
atoms and formation of overhangs and bulk defects is negligibly small. As a 
consequence the continuum model describing this processes must conserve the 
number of particles on the interface. The introduction of conservation laws into the 
growth equation forms new universality classes in surface phenomena. One of these 
classes is known as the Linear MBE equation (or the Mullins-Herring (MH) 
universality class [15]) and is described, in Fourier components, by the equation 



^ = -Kq 4 h q+ Tl q (t), (1) 

where h q is the Fourier component of the height measured relative to its spatial 
average, and t) (t) is the fluctuation of the rate of deposition, which is assumed to 
have a Gaussian distribution with zero mean and generally satisfies 

(^(t)^( t ')) = 2D q- 2p 5 q ,_ q ,5(t-f), (2) 

where D is a constant, p is a parameter that can be either positive or negative 
(actually in the case of conserved noise, which is of great interest in the MBE system, 
p = -1 , see re. [16]), and 8, is just the Kronecker symbol. 

However, just like in the non-conservative case, a non-linear extension was 
needed to describe the richness of the MBE processes. Various symmetry arguments, 
originally suggested by Villain [4, 20] as well as some physical arguments [1] indicate 
that the relevant MBE growth equation is 

A = -Kq 4 h q - -JL £q 2 (£ ■ m)5 q ,, +m h,h m + ^ , (3) 

dt V" e,m 

where g is the coupling constant, Q. is the volume of the system, to be taken 

eventually to infinity, and r\ is the noise term. This equation, which is by no means 

trivial, has been analyzed later using the Dynamical Renormalization Group method 
(DRG) [6]. The theoretical predictions for the critical exponents made, using this 
method, agreed quite well with results of numerical integration of eq. (3) as well as 
with results of simulations of discrete models belonging to the MBE universality class 
(see for example [1], [15] and references therein). Therefore, these results became 
widely accepted in the community of surface-growth physicists. 

However, in recent years some researchers raised again the question of the 
validity of the DRG theoretical predictions. One line of criticism was taken by 



DasSarma [15] who pointed out that the DRG results are derived from a leading order 
£ - expansion of a one-loop renormalization analysis, where 8 = 4 - d (d being the 
substrate dimension). He stressed the point that for the relevant dimensions discussed 
in the literature, i.e. d = 1 or 2 , the expansion parameter 8 = 3 or 2 is not small, 
therefore one may legitimately question the validity of the calculated exponents. 

A somewhat different line of criticism, however more radical and explicit was 
taken by Janssen [16] who was able to show that a two-loop calculation gives non 
trivial (although small) corrections to the critical exponents predicted by one-loop 
DRG calculation (more specifically oc two _ loop = oc one _ loop - 8 , where a is the roughness 

exponent and 8 is the small correction). By doing so he actually made a substantial 
contribution in refuting the underlying assumption that the coupling constant 
renormalizes trivially. This assumption was very essential to the one-loop DRG 
calculations done so far. Janssen was also able to show explicitly the reason for this 
discrepancy, which has to do with a mathematically ill defined generalization of so 
called the Galilean invariance (actually tilt invariance) of the KPZ equation [3] 
suggested by Sun, Guo & Grant [17]. As mentioned above, Janssen found that the 
correction to the scaling exponents was very small, and he suspected that the 
smallness of the correction was related to many, but incomplete, cancellations 
between diagrams as well as within internal momentum integrals. On that basis he 
speculated that a mode-coupling approach is a useful approach for this problem. 

In this paper I apply a method developed by Schwartz and Edwards [11-13] 
(also known as the Self-Consistent-Expansion (SCE) approach). This method has 
been previously applied to the KPZ equation. The method gained much credit in being 
able to give a sensible prediction for the KPZ critical exponents in the strong coupling 
phase, while DRG was not able to give any prediction for that phase for d > 1 (only 



the weak coupling solution was addressed). It is worth mentioning that this method is 
closely related to the mode-coupling approaches, in the sense that similar (but not 
identical) equations are obtained, while the underlying derivation is different. It is 
therefore, hoped that this paper will help to decide this unresolved situation, thus 
facing the challenge set up by Janssen [16]. I obtain the original results of the 
one-loop calculation [6, 8], thus corroborating these results, while avoiding the 
mathematical pathologies faced by the DRG method. This situation where DRG 
results obtained from different orders in 8 give different conceptual scenarios (i.e. 
trivial Vs. nontrivial renormalization of the coupling constant) calls for a resolution. 

Another remarkable advantage of the SCE method is the minor changes 
needed in order to generalize the result with uncorrelated noise to include noise 
correlated in space. The above implies a second important motivation for this paper, 
namely a demonstration of the robustness of the SCE method as well as its 
mathematical coherence and consistency. 

The SCE method is based on going over from the Fourier transform of the 
MBE equation in Langevin form to a Fokker-Planck form and constructing a 
self-consistent expansion of the distribution of the field concerned. The expansion is 
formulated in terms of (|> and oo q , where <|> is the two-point function in momentum 



space, defined by (|) q = (h q h_ q j , (the subscript S denotes steady state averaging), and 
co q is the characteristic frequency associated with h q , defined by 

<=j(h q (t)h_ q (o))dt/v 

o / 

I expect that for small enough q, (j) q and oo q are power laws in q, 

(b = Aq~ r and CO = Bq 1 * (4) 



[Since dynamic surface growth is a remarkably multidisciplinary field, there are 
almost as many notations as there are workers in the field. Therefore I give a brief 
translation of our notations to those most frequently used: 

|i = z, cc=(r-d)/2, and p = a/z = (r-d)/2ji. (5)] 

The method produces, to second order in this expansion, two nonlinear 
coupled integral equations in <j) and oo q , that can be solved exactly in the asymptotic 

small q limit to yield the required scaling exponents governing the steady state 
behavior and the time evolution. 

I begin with writing the Fokker- Planck form of the MBE equation (eq. (3)) 



— + > 

at v^ q 



D o q ^-+K q h q+ £ Mq , m hA 



g „2, 



P = 0, 



(6) 



where K q = Kq 4 , D 0q =q 2p and M qtm = -^q 2 (i. m )5 q ,, + „ 



A self-consistent expansion for such an equation was derived in the past (ref. 
[11-13]). The main idea is to write the Fokker-Planck equation 3P/3t = OP in the 
form 3P/3t = [O +0! +0 2 ]P , where O is to be considered zero order in some 
parameter X, O, is first order and 2 is second order. The evolution operator O is 



chosen to have a simple form O = -V 



D <ahT +0 ^ 



, where D q /w q = 4> q . 



Note that at present (|) and oo q are not known. I obtain next an equation for the 
two-point function. The expansion has the form (|) = <)) q + c \(j) p ,co p j, because the 
lowest order in the expansion already yields the unknown (|) . In the same way an 
expansion for oo q is also obtained in the form co q = oo q +d q {(|) p ,(JO p }. Now, the 



two-point function and the characteristic frequency are thus determined by the two 
coupled equations 

c q fe P ,« p }=0 and d q fo p ,G> p }=0. (7) 

Working to second order in the expansion, one gets the two coupled integral equations 

U 0q -^q% +I L 2 2- 2 2- = U ' W 



l,m « q + (X) t + (0„ 



W q + (0, + w„ 



W q + CO, + w„ 



K q -co q -2£M 



M, mq (^ m +M m , q ^ 



q^m 



= 0. 



C0,+CQ n 



(9) 



where in deriving the last equation I have used the Herring consistency equation [14]. 
In fact Herring's definition of oo q is one of many possibilities, each leading to a 

different consistency equation. But it can be shown, as previously done in [12], that 
this does not affect the exponents (universality). 

In the following I will treat the equations (8) and (9) for our specific problem 
of interest (i.e. use the specific form of K q and M qlm for MBE). These equations can 

be rewritten as 



D q- 2p -Kq> q + I 1 (q)(^ q +I 2 (q) = 
C0 q -Kq 4 +J(q) = 0, 



(10) 
(ID 



where the functions I, (q) , I 2 (q) and J (q) are given by 



i 1 (q) = 



2g 2 r d q 2 p-(q-^) 

■^— \d d £ L V L 

(2%) J co q + co,+co q _ 



q-l\ 2 £-q^ + £ 2 (q-iyq^_ t 



, (12) 



I 2 (q) : 



2g 2 



w 



yt. 



p.(q-i); 



C0 q +CQ,+C0 q _, 



-<M>q- 



%)=■ 



2g 2 r d q 2 p-(q-^) 

^-r\d d £ L V 

(2jc) j ©, + ©,_, 



q_ip.q^ + ^(q_i).q^ 



(13) 



(14) 



As previously stated (eq. (4)) I expect that for small enough q, (b and oo q are 
power laws in q, for small q (i.e. (b = Aq~ r and co q = Bq^). I am interested in 
eqs. (10) and (11) for small q's only. But, in order to achieve that one must consider 
the contribution of the large £ integration on the small q behavior of the whole 
integrals. So I break up the integrals I ; (q) and J(q) into the sum of two 

contributions I*(q), J > (q) and I*(q), J < (q) , corresponding to domains of £ 
integration, with high and low momentum respectively. I expand Ij' (q) and J* (q) 

for small q's and obtain the leading small-q behavior of the integrals, and after 
retaining only the leading terms, eqs. (10) and (11) reduce now to 

D q- 2p + A 2 q 4 - Kq> q + A iq 6 (b q + 1< (q)<b q + 1< (q) = (15) 

(O q -Kq 4 +J < (q) + A 3 q 6 =0. (16) 

At the mere price of renormalizing some constants in both equations, I am left 

with the integrals Ij" (q) , V 2 (q) and J < (q) that can be calculated explicitly for small 

q's since for small £ ' s the power law form for <b^ and (O t (or (b t and (H q _ t ) can be 

used (eq. (4)). In addition, the small-q dependence of each of the integrals naturally 
depends on the convergence of the integrals without cutoff. So, to leading order in q 

</ x </ x k 6 for d+4-r-(i>o 

I^(q),r(q)oc 4 ** (17) 

lKJ KJ lq d+8 ^ for d + 4-r-^t<0 

<t , fq 4 for d+4-2r-(l>0 

n(q)oc 4 ^ (18) 

2V J \ q ^-» for d + 4-2r-^t<0 
Going through the steps of a detailed analysis (like in the appendix of ref. 
[12]) I find that above d c =4+2p a weak-coupling solution, where the exponents 



r = 4+2p and p = 4 are obtained. The lower critical dimension of the non-linear 
MBE equation is thus 4 + 2p . 

A strong coupling solution can be obtained provided the T and p obey 
d + 4-T-p<0 and d + 4-2r-p< . In that case eqs. (15) and (16) take the form 

D q- 2 ' + A 2 q 4 - KAq 4r + AA.q^ + -^ r ^q d+ —F(r,p) = , (19) 

(2%) & 

ll 1 A 
Bq^ - Kq 4 + A,q 6 + ^L_^ q **"^G (r,p) = . (20) 

(271 ) o 

where F(T,p) is given by 

t-(e-t) 



F(r,p) = -fd d t — ^ — L Tie _ 1 1 2 ( t ■ e) t- r + 1 2 (e - 1 ) ■ a la - 1 r r 

t** + ie — t r + 1 L 

rt-(e-t)] 2 



+ 

(21) 



+ 

t^+le-tr + 1 



and G (r,p) is given by 



G(r,p) = -[d d t ^ ^ |e-t| 2 (t-e)r r + t 2 (e-t)-e|e-t| 



t^+le-tl 



(22) 



From the conditions given above for a strong-coupling solution I find 
d + 8-r-p<4. Therefore, the last term in eq. (20) is dominant over the second 
term. Two possibilities seem to arise now. Either the last term dominates the first term 
in eq. (20), which implies G(T,p) = 0, or these terms are proportional to the same 

power of q, which implies the scaling relation d + 8-r-2p = 0. The first possibility 
requires d + 8-r-2p<0, which is inconsistent with the whole idea of the 
expansion. The point is that higher order corrections have additional powers of 
qd+8-r-2n sQ t ^ at | n Qur case t ^ e j-gqyjj-ejnent d + 8-r-2p<0 means that higher 



order terms are more violent than lower order ones (for small q's ). Such a situation 
either implies inconsistency of the expansion, or calls for summing up the whole 
series in order to get a meaningful result. I assume that the expansion is consistent so 
that I am left with the second possibility, i.e. d + 8 -r-2|i =0 . It is interesting to 
mention here that Janssen's result is also consistent with the assumption that the 
requirement d + 8-r-2(i<0 cannot be fulfilled. 

As for eq. (19), I get d + S-2T-\x<4-T, meaning that the third term is 
negligible compared to the last. Here again, I am faced with two possibilities. Either 
the first term and the last term are proportional to the same power of q, resulting in 
d + 8-2r-(i = -2p , or the last term is dominant over the first term, in which case 

d + 8-2r-(i< -2p , and F(r,|j,) = 0. A careful numerical calculation shows that 
F(r,|i) ■$■ for d < 4 , so that the second possibility is ruled out. Therefore, I am left 

with the first possibility, namely with d + S-2T-\x = -2p . 

I am led to the conclusion that the strong coupling solution is 
r=(d + 4p + 8)/3 and |0,= (d-2p + 8)/3 . Taking into account the condition for a 

strong-coupling solution given above, I find that such a solution is valid only for 
d < d c = 4 + 2p . For d > d c the exponents describing the system are the same as those 
of the linear MBE equation, i.e. T = d and (I = 4 . It should be mentioned that The 
p- dependent strong coupling solutions exist as long as p>(d-4)/2 (for lower 
values of p we are actually always below the critical dimension, and the critical 
exponents are consequently always T = d and (1 = 4). The translation of these results 
to the frequently used notation reads 

z = j(d-2 P + 8)/3 d<4 + 2p 

[4 d>4+2p 

10 



(2p + 4-d)/3 d<4+2p 
a= v " . (24) 

[0 d > 4 + 2p 

The final conclusion is that the second order self-consistent expansion yields 
results that corroborate the results of one-loop DRG [6, 8]. As mentioned at the 
beginning of the paper, the mode coupling approach give similar equations to those of 
SCE to second order - although a different derivation and analysis (i.e. summation of 
the perturbation series while neglecting vertex renormalization Vs. a perturbation 
theory for the Fokker-Planck form). On that basis, I expect the same results from a 
mode-coupling approach when applied to this problem. 

To evaluate these results two facts have to be taken into account. The first is 
that the SCE approach does not rely on the symmetry argument [17] that its possible 
weakness was pointed out by Janssen. The second is that the SCE is known in other 
cases (i.e KPZ, see refs. [11-13]) to yield results that deviate from the results of the 
one-loop DRG and agree much better with simulations. This suggests that the 
one-loop DRG results may be exact for the case of MBE. On the other hand, attempts 
to verify Janssen's results via numerical simulations indeed found such corrections 
[18-19]. However, the corrections were systematically much larger than those 
predicted by Janssen himself. For example, for d = 1 the deviation 8 from the 
one-loop result differs from Janssen's result by an order of magnitude (correction of 
8=0.0025 Vs. 8 = 0.02). Since 8=0.0025 is closer to 8=0 than to 8 = 0.02, this 
suggests that the deviation found in the numerical simulations might be related to 
some other factors. Clearly further investigation is needed to decide this interesting 
issue. 
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